home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / CmplxQRP.cc < prev    next >
C/C++ Source or Header  |  1997-03-07  |  3KB  |  154 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #if defined (__GNUG__)
  24. #pragma implementation
  25. #endif
  26.  
  27. #ifdef HAVE_CONFIG_H
  28. #include <config.h>
  29. #endif
  30.  
  31. #include <cassert>
  32.  
  33. #include "CmplxQRP.h"
  34. #include "f77-fcn.h"
  35. #include "lo-error.h"
  36. #include "mx-inlines.cc"
  37.  
  38. extern "C"
  39. {
  40.   int F77_FCN (zgeqpf, ZGEQPF) (const int&, const int&, Complex*,
  41.                 const int&, int*, Complex*, Complex*,
  42.                 double*, int&);
  43.  
  44.   int F77_FCN (zungqr, ZUNGQR) (const int&, const int&, const int&,
  45.                 Complex*, const int&, Complex*,
  46.                 Complex*, const int&, int&);
  47. }
  48.  
  49. // It would be best to share some of this code with ComplexQR class...
  50.  
  51. ComplexQRP::ComplexQRP (const ComplexMatrix& a, QR::type qr_type)
  52.   : ComplexQR (), p ()
  53. {
  54.   init (a, qr_type);
  55. }
  56.  
  57. void
  58. ComplexQRP::init (const ComplexMatrix& a, QR::type qr_type)
  59. {
  60.   assert (qr_type != QR::raw);
  61.  
  62.   int m = a.rows ();
  63.   int n = a.cols ();
  64.  
  65.   if (m == 0 || n == 0)
  66.     {
  67.       (*current_liboctave_error_handler)
  68.     ("ComplexQR must have non-empty matrix");
  69.       return;
  70.     }
  71.  
  72.   Array<Complex> tau (m < n ? m : n);
  73.   Complex *ptau = tau.fortran_vec ();
  74.  
  75.   int lwork = 3*n > 32*m ? 3*n : 32*m;
  76.   Array<Complex> work (lwork);
  77.   Complex *pwork = work.fortran_vec ();
  78.  
  79.   int info = 0;
  80.  
  81.   ComplexMatrix A_fact = a;
  82.   if (m > n)
  83.     A_fact.resize (m, m, 0.0);
  84.  
  85.   Complex *tmp_data = A_fact.fortran_vec ();
  86.  
  87.   Array<double> rwork (2*n);
  88.   double *prwork = rwork.fortran_vec ();
  89.  
  90.   Array<int> jpvt (n, 0);
  91.   int *pjpvt = jpvt.fortran_vec ();
  92.  
  93.   // Code to enforce a certain permutation could go here...
  94.  
  95.   F77_XFCN (zgeqpf, ZGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork,
  96.                  prwork, info));
  97.  
  98.   if (f77_exception_encountered)
  99.     (*current_liboctave_error_handler) ("unrecoverable error in zgeqpf");
  100.   else
  101.     {
  102.       // Form Permutation matrix (if economy is requested, return the
  103.       // indices only!)
  104.  
  105.       if (qr_type == QR::economy)
  106.     {
  107.       p.resize (1, n, 0.0);
  108.       for (int j = 0; j < n; j++)
  109.         p.elem (0, j) = jpvt.elem (j);
  110.     }
  111.       else
  112.     {
  113.       p.resize (n, n, 0.0);
  114.       for (int j = 0; j < n; j++)
  115.         p.elem (jpvt.elem (j) - 1, j) = 1.0;
  116.     }
  117.  
  118.       if (qr_type == QR::economy && m > n)
  119.     r.resize (n, n, 0.0);
  120.       else
  121.     r.resize (m, n, 0.0);
  122.  
  123.       int min_mn = m < n ? m : n;
  124.  
  125.       for (int j = 0; j < n; j++)
  126.     {
  127.       int limit = j < min_mn-1 ? j : min_mn-1;
  128.       for (int i = 0; i <= limit; i++)
  129.         r.elem (i, j) = A_fact.elem (i, j);
  130.     }
  131.  
  132.       int n2 = m;
  133.       if (qr_type == QR::economy)
  134.     n2 = min_mn;
  135.  
  136.       F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
  137.                  pwork, lwork, info));
  138.  
  139.       if (f77_exception_encountered)
  140.     (*current_liboctave_error_handler) ("unrecoverable error in zungqr");
  141.       else
  142.     {
  143.       q = A_fact;
  144.       q.resize (m, n2);
  145.     }
  146.     }
  147. }
  148.  
  149. /*
  150. ;;; Local Variables: ***
  151. ;;; mode: C++ ***
  152. ;;; End: ***
  153. */
  154.